ANALYSIS OF THE MONTE-CARLO ERROR IN A HYBRID 
SEMI-LAGRANGIAN SCHEME 



CHARLES-EDOUARD BREHIER AND ERWAN FAOU 

Abstract. We consider Monte-Carlo discretizations of partial differential equations based on a 
combination of semi-lagrangian schemes and probabilistic representations of the solutions. We 
study the Monte-Carlo error in a simple case, and show that under an anti-CFL condition on the 
time-step St and on the mesh size 5x and for A'' - the number of realizations - reasonably large, we 
control this error by a term of order 0{yj5t/N). We also provide some numerical experiments to 
confirm the error estimate, and to expose some examples of equations which can be treated by the 
numerical method. 



1. Introduction 

The goal of this paper is to analyze and give some error estimates for numerical schemes combining 
the principles of semi-lagrangian and Monte-Carlo methods for some partial differential equations. 
Let us describe the method in a very general case by considering first a linear transport equation 
of the form 

dtu{t,x) = f{x) ■\7u{t,x), xGM.'^, u{0,x) = uo{x), 

where uq is a given function. Under some regularity assumptions and existence of the flow associ- 
ated with the vector field f{x) in M'^, the solution of this equation is given by the characteristics 
representation u{t, x) = u{0, ipt{x)), where ^t{x) is the flow associated with the ordinary differential 
equation y = f{y) in M'^. In this context, semi-Lagrangian schemes can be described as follows. 
Let us consider a set of grid nodes Xj, j ^ K in M*^ [K = N or a finite set) and an interpolant 
operator I mapping vectors of values at the nodes, (uj) G to a function (Iu){x) defined over 
the whole domain. In this paper we will consider the case where Xj = j{6x), j G Z'^, 6x is the space 
mesh size, and X a standard linear interpolation operator. Given approximations Uj of the exact 
solution u{tn,Xj) at times t„ = n{6t) and points xj, the previous formula gives an approximation 
scheme for u^~^^ obtained by solving the ordinary differential equation y = f{y) between t„ and 
tn+i- ^j"*"^ = {^u"'){^st{xj)), where is the numerical flow associated with a time integrator. 

These methods are particularly interesting when the vector field f(x; u) depends on the solution 
u making the transport equation nonlinear, see for instance |12| HI [11] and the references therein. 
This is the case when an advection term is present for instance, or for Vlasov equations (see for 
instance [2]). In these situations, standard semi-lagrangian schemes are based on solving equations 
of the form 

dtu{t, x) = /(x; u") • Vu{t, x) , 

between tn and tn+i, where denotes the solution at time t„. In other words the vector field is 
frozen in (in the language of geometric numerical integration, it is Crouch and Grossman method, 
see ^j). If moreover the vector field f{x;u) possesses some geometric structure for all functions u, 
the numerical integrator can be chosen to preserve this structure (for example symplectic integrator 
in the Vlasov case). 
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In many situations, a diffusion term is present, and the equation can be written (in the Hnear 
case) 

(1.1) dtu{t,x) = ^a{x)a{x)'^Au + f{x) ■Vu{t,x), xGR'^, u{0,x) = uo{x), 

where ct(x) is a d x n matrix. In this case the solution admits the characteristic representation 

u{t,x) = Etio(Xf), 

where is the stochastic process associated with the stochastic differential equation 

(1.2) dXf = /(Xf)dt + a(Xf)d5t, XS = x, 

where {Bt)t>o is a standard n-dimensional Brownian Motion. 

In general, the law of the random variable Xf is not explicitly known, and we are not able 
to compute the expectation. The classical approximation procedures for such problems are Monte- 
Carlo methods: if we assume that we are able to compute N independent realizations (-'^f ''")i<m<A'' 
of the law of , we can approach u{t, x) with 



1 ^ 

m=l 

In general, the variance of the random variables uo(-Yf '*") is of size t and the law of large numbers 
ensures that the statistic error made is typically of order 0{y^T/N) for an integration over the 
interval [0,T]. To this error must be added the error in the approximation of the process X^ 
by numerical schemes of Euler type for instance. This is error is of order 0{6t), see for instance 
[U [ini [l3] and the reference therein for analysis of the weak-error in the numerical approximation 
of stochastic differential equations. If a global knowledge of the solution is required, the above 
operation must be repeated for different values of Xj on the grid. 

The numerical method we study in this paper is based on the Markov property of the associated 
stochastic processes: we have for any Xj on the spacial grid and locally in time 

(1.4) n(t„+i,x,) =Eu(t„,Xj/), 

which is the formula we aim at discretizing. Using the Euler method to compute a numerical 
approximation of Xg^ , we end up with the following numerical scheme 

1 ^ 

(1.5) -■=mY1 (^^") (^i + + ^^AA-'-'^), 



N 



m= 



where the random variables {M^'"^'^)i<m<M are independent standard Gaussian random variables. 
Note that the main difference between the standard Mont-Carlo method is that the average is 
computed at every time step. 

Such numerical method were already introduced in [3j. The principle of using (random) charac- 



teristic curves in (1.4) over a time interval of size 6t and to use an interpolation procedure to get 
functions defined on the whole domain fits in the semi-lagrangian framework. The addition of the 
Monte-Carlo approximation then justifies the use of the hybrid terminology. 

As in the deterministic case described above, it is clear the method can be adapted in situations 
where the drift term f{x) and the noise term (t{x) depend on the solution u. We will present in the 
end of the paper some numerical experiment in such nonlinear situations. 

Another remark is that different kind of boundary conditions can be considered. If the presen- 
tation made above was concerning the situation where the equation is set on W^, representation 



formulae such as (1.4) hold true in the case of periodic boundary conditions, Dirichlet or Neumann 
condition on bounded domains. Again, we give some numerical examples in the end of the paper. 
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But the main aim of this paper is to perform the numerical analysis of the scheme (1.5) in the 
simplest situation, that is / = 0, in dimension d = 1, with crix) = 1 and periodic intial conditions 
such that the transport equation is set on a domain = (0, 1) C M with periodic boundary 



conditions. Note that using splitting strategy, this is the term that is new in (1.5) in comparison 
with standard semi-Lagrangian methods. 



In essence, the result stated in the next section shows that in this situation, the algorithm (1.5) 
approximates the exact solution up to an error that is of the order 0{6x + ^J5t/N). The first term 
comes from the interpolation, and the second show a variance reduction phenomenon. To obtain this 
bound, we require an anti-CLF condition as usual for semi-Lagrangian methods. We also assume 
that is sufficiently large in some relatively weak sense (see a precise statement below). 

The paper is organized as follows. In Section [2] we present the method, introduce various nota- 
tions and state our main result. In Section 3 we give some properties of random matrices arising 
in a natural way in the definition of the scheme ( |1.5[ ), and that are needed in the proof of our main 
estimate, given in full details Section 4. Possible extensions of the method, with for instance Dirich- 
let boundary conditions, or more complicated PDEs, are evoked in Section [5} together with a few 
numerical results. In particular, we present simulations for the two-dimensional Burgers equation. 

2. Setting and main result 

We consider the linear heat equation on (0, 1), with a smooth initial condition and with periodic 
boundary conditions; more precisely, we want to approximate the unique periodic solution of the 
following partial differential equation: 

^ ' dt 2 9x2' 

such that ti(0, .) = -uq ^ ~^ is a smooth periodic function of period 1, and such that for any 
t > the function u{t, .) is periodic of period 1. For periodic functions, we denote by and 
the usual function spaces associated with the respective norm and semi-norm 

\\f\\l2= flfix^dx, and \f\l,= [\f{x)\'dx. 
Jo Jo 

2.1. Interpolation operator. For a given integer Ms > 1, we discretize the space interval (0,1) 
with the introduction of nodes xj = j6x for j €z S := {0, . . . , Ms — 1}, with the condition xms — 
MsSx = 1. We set V = {{uj)jes} C M^^ the set of discrete functions defined on the discrete points 
of the grid. 

We use linear interpolation to reconstruct functions on the whole interval from values at the 
nodes. We define an appropriate basis made of periodic and piecewise linear functions for k £ S = 
{0, ...,M5-1}. We set for [xfc - 1/2, -M/2], 

<Pk{x) = ^ ^^ ), where (i>{x) = < ^ j^j ?^ ^ ^ 
ox I 1 — \x\ if |x| < 1, 

and extend the function (pk by periodicity on (0, 1). Hence (pk is the unique piecewise linear periodic 
function and satisfying (pk{xj) = 5kj the Kronecker symbol. Note that we have "^k^s ^k{x) = 1 for 
all x E (0,1). 

We define the following projection and interpolation operators 

^ (h^ V ^^^^ 

\f ^ {f{xj))jes \u = {uk)kes ^ T.k=(^^ ^k4>k{x), 

where denote the Sobolev space of periodic function on (0, 1). Clearly, P o X is the identity 
on V; nevertheless the distance between the identity and the composition of the operators X oV 
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depends on the functional spaces and on the norms. Below, we give the estimates that are useful 
in our setting. We just notice that Zo P(l) = 1, as '^j. (t>k{x) = 1 for all x. 

2.2. Discrete norms. For any u = {uj)j^s we define the discrete norm and semi-norm as 

II ii2 c- 2 1 I i2 c- (^7 + 1 ~ ^7')^ 

\H\p=Sx}_^Uj and \u\^,=Sx}^ Sr^ ' 
jes jes 

where we use the extension by periodicity of the sequence (uj) for the definition of the semi-norm: 

we thus have ums = ''^o- We also define a norm with = (||^t||£2 + I'^^l^i)"^^^- 

With these notations, we have the following approximation results: 

Proposition 2.1. There exists a constant c > such that for any mesh size 6x = 1/Ms, and any 
sequence u = (uj) eV we have: 

\u\fii = \Iu\jji, and ||tt||^2 = 112^^*11^2 + cSx^\u\f^i. 

Moreover, for any function f £ we have 

11/ - {XoV)f\\l. < c5x' (l/l^. + |(XoP)/|^,) . 

Proof . The first equality follows from a direct computation. The second one is proved expanding 
in the L? scalar product Tu = Uk4>k, and rewriting the sums in order to make the semi-norm 
appear: we have 

kes k,ees kes kes 

where we define um^ = uq and u-i = ums-1j that is we extend u E V hj periodicity. We also used 
the fact that for all k, 

{'f>k,4>e)L^=0^^^{k-'i-,k,k + l}, {(j)k,(l>k)L^ = and {<j)kAk-i)L^ = -^- 

Now, the equality contains HuH^i which appears with natural integration by parts - using peri- 
odicity: 

5x 
T 



1^*11^2 - II^^IIl2 = XI ~ ""fe-l) + "^kiUk - Uk+l)) 



keS 

= Y XI i'^k+iiuk+i - Uk) + Uk{uk - Uk+l)) 

fce5 

fce5 



To prove the last estimate, we have for any j € S = {0, 1, . . . Ms — 1} and for any x G [xj,Xj-^-l] 

\f{x) - (I o vf){x)f < 2(£ f'mf + 2(£ li^i+A^lM^tf 



< 2{x - Xj) j \f'it)\^dt + 2(5x^ 

J X-j 



\f{xj+i) - f{Xj)\ 
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Sx"^ ' 

using the Cauchy-Schwarz inequality. Now we integrate over [xj, Xj+i], and then it remains to take 
the sum over j G -S of the following quantities: 

(2.2) \f{x) - {XoVf){x)\^Ax < |/(x)|2dx + 25x2^^fell:^^^M^^^. 

Jxj Jxj ^X 
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The first term of the right-hand side is controlled with , while the second term involves I'P/l^i = 
\iIoV)f\l,. □ 

2.3. Definition of the numerical method. We consider a final time T > 0, and an integer Mt, 
such that we divide the interval [0, T] into Mt intervals of size 5t := We are thus interested 
in approximating the solution u{t,x) of (2.1) at times tn = nSt and nodes xj for n < Mt- In the 
simple case of Equation (2.1 ), for which we have the representation formula u{t, x) = Kuo{x + Bt), 
the numerical scheme (1.5) is the application n" i— t- n""*"^ from V to itself written 

(2.3) u]+' = ^^2(11 ^^^^'^i + ^AA"'-'^)] , 

m=l VfceS / 

where the random variables jV"'™'.?, indexed by < n < Mt — 1, 1 < m < N and j £ S 
are independent standard normal variables. More precisely, to avoid an error term due to the 
approximation of Brownian Motion at discrete times, we require that 

(2.4) V5tAf-'-- = B^-^l,^ - 5^7/-) 

for some independent Brownian Motions (S(™'J)) for 1 < m < iV and < i < Ms - 1. 

We start with an initial condition = {u^ = uo(xfc)), which contains the values of the initial 
condition at the nodes. To obtain simple expressions with products of matrices, we consider that 
vectors like are column vectors. 

We then define the important auxiliary sequence v'^ £ V satisfying the following relations: 

N 

n+1 _ _|_ 
^ N ^ , 

(2.5) ™=i ^kes 



with the initial condition = . Indeed, for any < n < Mt the vector is the expected value 
- defined component- wise - of the random vector u". 

2.4. Main result. With the previous notations, we have the following error estimate: 

Theorem 2.2. Assume that the initial condition uq is of class C^. For any p £ N and any final 
time T > 0, there exists a constant Cp > 0, such that for any 5t > 0, 6x > and N £N* we have 

5x'^ 

(2.6) sup \u{tn, Xj) — v^\ < C —-— sup |^io(^)l 

j&S Ot ^.g[o,i] 

and 

(2.7) E||." - v-\\% < C>°|^(1 + ^) (1 + I + g(l + I log(^OI))' (I + ^) • 

The control of the first part of the error is rather classical, while the estimate on the Monte-Carlo 
error given by (2.7) is more original and requires more attention in its analysis and in its proof. 

First, we observe that the estimate is only interesting if a condition of anti-CFL type is satisfied: 
for some constant c > we require 

Sx I 

— max(l, \/log((5t)) < c. 

We then identify in (2.7) a leading term of size which corresponds to the statistical error in a 
Monte-Carlo method for random variables of variance bt^ and a remaining term, which goes to 
with arbitrary order of convergence with respect to the number of realizations N . This second term 
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is obtained via a bootstrap argument. Indeed it is easy to get the classical estimate with p = 0. The 
core of the proof is contained in the recursion which allows to increase the order from p to p + 1; it 
heavily relies on the spatial structure of the noise and on the choice of the ^^-norm. 



Thanks to (2.7) when p = 1, we see that interpreting Theorem 2.2 as a reduction of variance to 



a size 6t is valid: we bound the error with 

6t 1 3 
iV ^ iV2 - 1" + 2iV2' 

which can be compared with a classical Monte-Carlo bound with the variance 6t: we have for any 
sample (Yi, . . . , l^v) of a random variable Y 



N - 2 2iV2' 



The control of the Monte-Carlo error in Theorem 2.2 relies on several arguments. Firstly, the first 
factor corresponds to the accumulation of the variances appearing at each time step - where two 
sources of error are identified: the random variables involve a stochastic diffusion process evaluated 
at time 5t, and an error is introduced by the interpolation procedure. To obtain another factor, we 
observe that the independence of the random variables appearing for different nodes implies that 



only diagonal entries of some matrices appear - see (4.12). However, this independence property 
also complicates the proof: the solutions are badly controlled with respect to the semi-norm. 
We then propose a decomposition of the error where the number of realizations N appears in the 
variance with the different orders 1 and 2: the first part is controlled by 5t and 5x, while the second 
one is only bounded. We finally use recursively this decomposition in order to improve the estimate, 
with a bootstrap argument. 



3. Random matrices 



The definition of the numerical scheme (2.3) can be rewritten with matrix notations: for column 
vectors of size Ms such that {u^)j = Wj, we see that 

(3.1) = P^'^^u", 
where the entries of square matrix satisfy for any 1 < j,k < Ms 

1 ^ 

(3.2) Pjj = - ^ M^j + ^/5^AA"'-'^■). 

m=l 

Moreover we decompose these matrices into independent parts: for 1 < m < 

1 ^ 

(3.3) " iV ^ 

m=l 

with the entries (P("'™))j,fc = M^j + V&tAf"'"^'^). 

We observe that the matrices are independent; in each one, the rows are independent, 

however in a row indexed by j two different entries are never independent, since they depend on 
the same random variable J\f^'"^'^-^ moreover, the sum of coefficients in a row is 1. 

All matrices pC"'™-) have the same law; we define a matrix Q = EP^'^'™') = EP("), by taking the 
expectations of each entry: for any j,k £ S 



(3.4) Q,, = E[<Pk{xj + V5tAA"'"^'^)]. 
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The right-hand side above does not depend on n, m since we take expectation. It only depends 
on j through the position xj, not through the random variable J\f^'"^'i , With these notations, the 



vectors f " satisfy the relation 7;"+^ = Qv"^ - see (2.5) - and we have for any n > 

n-l 

(3.5) ^" = n = • • • and = Q"u°. 

1=0 

We only present a few basic properties of the matrices p("''"^\ and Q. First, we show that 
they are stochastic matrices. Second, we control their behavior with respect to the discrete norms 
and semi-norms. In order to prove the convergence result, we need other more technical properties 
which are developed during the proof. 

Proposition 3.1. For any < n < — 1 and for any \ < m < N , almost surely pC"'™-) is 
a stochastic matrix: for any indices j,k E S we have Pj^'"^^ > 0, and for any j £ S we have 

\^ p(n,m,) _ I 

For any < n < Mt — 1, -P'-"-' is also a random stochastic matrix. 

The matrix Q is stochastic and symmetric - and therefore is bistochastic. 

Proof . The stochasticity of the random matrices is a simple consequence of their definition 

(3.2 ) and of the relations (pkix) > and X^^-g^ (t^kl{x) = 1. Since P^"^ is a convex sum of the pC"'"*), 



the property for those matrices also holds. 

Finally, by taking expectation Q is obviously stochastic; symmetry is a consequence of (3.4), and 
of the property 4>k+i{x) = (pkix — 5x): 



= E[Mxk - Xj - ^/5^AA"'"•J■)] = nMxk - Xj + V^AA"''™''^)] = Qkj, 

since 0o is an even function, and since the law of J\f"'"^'^ is symmetric and does not depend on j. 
However this symmetry property is not satisfied by the P-matrices, because the trajectories of these 
random variables are different when j changes. □ 

Thanks to the chain of equalities in the proof above, we see that Qj ^ only depends on k — j, but 
we observe that no similar property holds for the matrices p("-'™). 

We now focus on the behavior of the matrices with respect to the £^-norm. The following 
proposition is a discrete counterpart of the decreasing of the L^-norm of solutions of the heat 
equation. 

Proposition 3.2. For any < n < Mt — 1 and for any 1 < m < N , and for any u gV we have 



E||p("'™)u||22 < \\u\\% and E||p(")n||22 < ||n| 



2 



3.1 



Proof . According to the definitions above (3.1) and (3.3), we have for any index j {P^"'''^^u)j 
Sfee5 ^jfc'™^^^' Thanks to the previous Proposition 



we use the Jensen inequality to get 



E||p("''")u||22 = 5x Ve|(p("'™)u) 



<6xY,Y.^p^r^\u,\^^ 

jes kes 




now we use the properties of the matrix Q - it is a bistochastic matrix according to Proposition 3.1 
- to conclude the proof, since "^j^s Qj,^ = 1- The extension to the matrices P*-"-* is straightforward. 

□ 

7 



The matrix Q satisfies the same decreasing property in the £^-norm; moreover we easily obtain a 
bound relative to the h}-sevai norm: 



Proposition 3.3. For any u ^ V, we have ||Q'u||^2 < ||u||^2 and \Qu\fii < \u\ 



Proof . The proof of the first inequality is similar to the previous situation for the random 
matrices. To get the second one, it suffices to define a sequence u such that for any < j < Ms — 1 
we have uj = "^"^^ - with the convention ums = uq. Then thanks to the properties of Q we have 
Qu = Qu: for any j ^ S 



{6x)Quj = {Qu)j+i - {Qu)j = ^Qj+i,kUk - ^Qj,kUk 

kes kes 

= ^ Qj,k-iUk - ^ Qj,kUk = ^ Qj,k{uk+i - Uk) = 6x{Qi 



kes kes kes 

using a translation of indices with periodic conditions, and the equality Qj^i^k = Qj,k-i ^ explained 
above. As a consequence, we have = ||(5m||£2 = ||(5m||£2 < ||n||£2 = |u|/ji. □ 

It is worth noting that the previous argument can not be used to control K\P^'^'"^^u\fii: for a 
matrix P = p(^'^\ the corresponding quantity Pu can not be easily expressed with u. Indeed, 
given a deterministic u, then {P^'^'"^^u)j and (p("'''"^ii)j_|_i are independent random variables - since 
they are defined respectively with J\f("'"^'^) and A^(">™''i+i). The only result that can be proved is 



Proposition 3.4 below. However, its only role in the sequel is to explain why we can not obtain 



directly a good error bound; as a consequence, we do not give its proof. 

Proposition 3.4. There exists a constant C, such that for any discretization parameters N > 1, 
6t = and 6x = -pj, we have for any vector u €z V 



(3.6) IE|P(%|^<(1 + C^y^)|n|^.. 



Due to independence of matrices involved at different steps of the scheme, the previous inequalities 
can be used in chain. 

We thus observe that the matrices P^^^ and Q are quite different, even if Q = EP'^^l On the 
one hand, the matrix Q is symmetric, and therefore respects the structure of the heat equation - 
the Laplace operator is also symmetric with respect to the L^-scalar product. On the other hand, 
the structure of the noise destroys this symmetry for matrices P^''\ while it introduces many other 
properties due to independence - in some sense noise is white in space and implies first that solutions 
are not regular, but that on the average a better estimate can be obtained. 



4. Proof of Theorem 12.21 



We begin with a detailed proof of (2.7). A proof of the other part of the error (2.6) is given in 
Section 4.4 below. Easy computations give the following expression for the part corresponding to 
the Monte-Carlo error: for any < n < Mt 



Ms-l Ms-1 

5x Var(u5') = ^ E| 
i=o i=o 



where the superscript * denotes transposition of matrices. 



Since the vectors u" and satisfy (3.5), with the same deterministic initial condition u^, we 
have 

EWu^'-v^l =E||(P("-1)...P(°) -Q")n°||22 

= (5x('U°)*E (^(p("-i) . . . - Q")*(p("-i) . . . p(0) - 
= (5x('u°)*E ((p(°))* . . . (p("-i))*p(«-i) . . . p(o) _ (g")*g") 



where the last inequality is a consequence of the relation EP*^'^) = Q and of the independence of the 
matrices P^*^^. 

Therefore we need to study the matrix = E ((p(°))* . . . (p{'^-i))*p{"-i) . . . p(o) _ (Q")*Q") 
given by the expression above, such that 

IE||n"-f"||^2 =(5x(n°)*5„,u°. 

4.1. Decompositions of the error. We propose two decompositions of Sn into sums of n terms, 
involving products of matrices P^^\ of Q and of the difference between two matrices P^^^^ and Q, 
which corresponds to a one-step error: 

n-l 

(4.1) P("-1) . . . P(°) - Q" = ^ p("-i) . . . p(fe+i) (pC^) _ Q)Qfc, 

fc=0 

and 

n-l 

(4.2) p{«-l) . . . p(0) _Qn QU-^l-k (^p{k) _ Q)p(fc-l) , _ _ p{0)_ 

fc=0 

These decompositions lead to the following expressions for Sn - where we use the independence 
of the matrices P^'^) for different values of k: 

n-l 

fe=0 
n-l 

fc=0 

Therefore we obtain the following expressions for the error: 
E\\u'' - = Sxiu'^y SnU° 

= y E||p("-i) . . .pC^+DrpW - g)g'=«0||2 

(4.3) to 

n-l 

= ^E||Q"-i-'=(pW - Q)P('=-1) . . . P(°)n0||22. 

fc=0 

Before we show how each decomposition is used to obtain a convergence result, we focus on the 
variance induced by one step of the scheme. In fact, only the second one gives the improved estimate 
of Theorem |2.2[ Nevertheless, we also get a useful error bound thanks to the first one. 
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5„ = E Y.iQ'^T [P^''^ - Q) {P^^^^'^T . ■ ■ (p("-i))*p("-i) . . . p(^'+i) (^p(fc) - Qj q'^ 

fe=0 
n-l 

= E^(p(°))* . . . (p(^-i))* ('p('=) - qV (g"-i-'=)*Q"-i-'=('pW - q)p(^"1) . . . p(o) 



4.2. One-step variance. In the previous Section, we have introduced decompositions of the error, 
and we observed that we need a bound on the error made after each time-step. The following 
Proposition states that the variance after one step of the scheme is of size 6t if we consider the 
£^ norm, and that a residual term of size 6x'^ appears due to the interpolation procedure. If we 



consider independent realizations. Corollary 4.2 below states that the variance is divided by 
if we look at the full matrix of the scheme. 

Proposit ion 4.1. There exists a, constant C, such that for any discretization parameters 6t — 
and 6x = and for any 1 < m < N and < n < Mt — 1, we have for any vector u G 

(4.4) E||(p("''")-Q)^x||22 <C{5t + 5x')\u\l,. 

Corollary 4.2. For any < n < Mt — 1 and for any vector u G M*^^, we have 



1 • 



The proof of the corollary is straightforward, since = ^ with independent and 



identically distributed matrices However, the proof of Proposition 4.1 is very technical. 

One difficulty of the proof is the dependence of the noise on the position j: for different indices 
ji and j2, the random variables (p("''"^)ii)j^ and {P^"'''^^u)j2 are independent. To deal with this 
problem, for each j we introduce an appropriate auxiliary function and we analyze the error on 
each interval separately. We also need to take care of some regularity properties of the 

functions - they are functions, piecewise linear, but they are not in general of class - in order 
to obtain bounds involving the and semi-norms. 



Proof of Proposition 4.1 To simplify the notations, we assume that n = and that m = Iso 



that we only work with one matrix P with entries 

Pj,k = 4>k{xj + -BjJ, 

where the are independent Brownian Motions. 

We define the following auxiliary periodic functions: for any x G 

(4.5) V{x) =EIu{x + B^^), 
and for any index < j < Ms — 1 

(4.6) U^^\x) =Iu{x + Bi^). 



We observe that since we take expectation in (4.5) the index j plays no role there. Moreover we 
have the following relations for any j G S: 

V{xj) = {Qu)j and u'^^Xxj) = {Pu)j, but U^j\xj+i) ^ {Pu)j+i. 

The last relation is the reason why we need to introduce different auxiliary functions U^^^ for each 
index j. 

We finally introduce the following function depending on two variables: for any < t < 5t and 

x G M, 

(4.7) Vit,x) =EIu{x + Bt), 

for some standard Brownian Motion B. This function is solution of the backward Kolmogorov 
equation associated with the Brownian Motion, with the initial condition V(0, .) = lu, and for 
t > 

dtV = \dl,V. 

10 



Moreover we have V{5t, .) = V. 

We have the following expression for the mean-square error, integrated over an interval 
for any index j E S 



(4.8) 



E\U^^\x) -V{x)\'^dx 



St px 



Cj + l 



E\d^V{St - s,x + Bi)\'^dxds. 



The proof of this identity is as follows. First, thanks to smoothing properties of the heat semi-group, 
for any t > the function V{t, .) is smooth. Using Ito formula, with the Brownian Motion B^^^ 
corresponding to the function U^^\ 

dV{6t -s,x + Bi) = d^V{5t -s,x + Bi)dBi, 

for < s < 5t — € and for any e € (0, 5t), and the isometry property implies 



E\V{6t,x)-V{e,x + Bi ^'^ 



St- 



\d^Vi6t- s,x + Bi)\'^ds. 



We integrate over x G [xj,Xj^i], and we then pass to the limit e — )• 0, since V(0, .) = lu is a 
piecewise linear function. Moreover, we use the identity V{6t, .) = V. We observe that in the 
right-hand side of the last equality we take expectation, so that we replace B^ with the Brownian 
Motion B, which does not depend on j. 

Summing over indices j £ S, we then get, thanks to an affine change of variables y = x + Bg 



St rl 



V / ^U'^^\x)-V{x)\^dx= / / 

,^cJx^ Jo Jo 



Jo 

St f\ 



E|axV((5t-s,x + S,)pdxds 

■St 



pot pi pot 

= / \d^V{6t-s,x)\'^dxds= / \Vi5t- s,.)\jjids 
Jo Jo Jo 

rSt 

< / |V(0, .)|^ids = 6t\Iu\l, = dt\u\l,. 
Jo 



The inequality (|4.4|) is then a consequence of the two following estimates: first. 



(4.9) 



and second we show that 



Xj + l 



E 



dx < C6x'^\u\f^i , 



(4.10) 



E\\Pu-Qu\\% [ E\IoV{U^^^ -V){x)\'^dx <C6x'^\u\li. 

c' J Xi 



To get (4.9), we use the inequality (2.2) on each interval for a fixed realization of i?^^: 



Xj + l 



|C/(^')(a;) - y(x) - Xop([/W - V){x)\'^dx < Cbx'^ 



Xj+l 



-|- C5x6x 



2 |[^(^-)(x,+i) - V{x,+^)] - [U(^\x,) - V{x,)]\^ 



Sx"^ 



Taking the sum over indices j £ S and expectation, we see that 

E|a^.([/(^) - V){x)\'^dx < 2^ E\d^{Iu){x + Bi)\'^dx + ^ H''' \d^V{x)\'^dx 



je5 -^3 jeS 

1 2 I \|2 \ ^ /il'7'-,,|2 /1 1^,1 2 



< 2i\Iu\i,i + \V{5t, .)|^0 < M^Hi = Mu\ii, 

since V = V{5t,.). Indeed, taking expectation allows to consider a single Brownian Motion B, 
without j-dependence. 
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We now decompose the remaining term as follows: 



+ 2 



6x'^ 

\V{x,+,)-V{x,)\^ ^ 
5x'^ 



With the second part, using Proposition 3.3 we see that 
\Vix,+,)-V{x,)\' 

ie5 



Sx^- 



Sx"^ 



^ bx^ 



\Qu\\^ < \u\li. 



To treat the first part, we make the fundamental observation that for a fixed j E S, the same noise 
process is used to compute all values U^^\x) when x varies. As a consequence, we can use a 
pathwise, almost sure version of the argument leading to the proof of Proposition 3.3 which concerns 
the behavior of Q with respect to the semi norm. 

U^^Hxj+i) - U^^\x,) = Uk[<Pk{xj+i + Bl) - M^j + Bi)] 

= M<Pk-iixj + Bl^) - (pkixj + Bl^)] 
k&S 

= - Uk]4'k{Xj + 

kes 

using the relation <j)k+i{x) = (pkix — Sx) and an integration by parts. 

Now summing over indices j G 5 and using the Jensen inequality - thanks to Proposition 3.1 
we obtain 



|2 



dx"^ 



<6xY 

k&S 



\Uk+l - Uk\ 

dx"^ 



Having proved (4.9), we now focus on (4.10). We have, since 

_ v)j = [{P - Q) U\j 

Y \I oV{U^^'^ -V){x)\^dx - SxYlKP - Q>]j\'^\ 

|(c/0-)-y)(x,+i)-(t/0-)-y)(x,)P 



C5x^5x Y 



dx"^ 



It remains to take expectation and to conclude like for (4.9). 



□ 



4.3. Proof of Theorem 2.2, As we have explained in the introduction, we consider that 5x is 
controlled by 5t thanks to a anti-CFL condition. Roughly, from Proposition |4.1| we thus see that 
the variance obtained after one step of the scheme is of size and that the error depends on the 

and 



3.3 



3.4 



solution through the h} semi-norm. Moreover, from Propositions 
behaviors of the matrices Q and P*^") with respect to this semi-norm are quite different. 
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we remark that the 



Using the first decomposition of tlie error in (4.3), we use in chain the bounds given above in 



Propositions 3.2 3.3 and 4.1 and Corollary 4.2 

n-1 



n-1 



k=0 



k=0 



k=0 



If the continuous problem is initialized with the function uq, which is periodic and of class C^, 
then = Vuq satisfies < sup^gjg,!] l''^o(^)l- Moreover we assume that an anti-CFL condition 

is satisfied, so that the term 5x^ /5t is bounded. As a consequence, we find a classical Monte-Carlo 
estimate, where the error does not decrease when 6t goes to and is only controlled with the number 
of realizations: 



(4.11) 



E|K-»"l&<ci±^|.»li.. 



In fact, (4.11 ) shows that the variances obtained at each time step can be summed to obtain some 
control of the variance at the final time. To get an improved bound, we thus need other arguments. 

The main observation is that using independence of rows in the P-matrices, we only need to focus 
on diagonal terms 



sup((QTQ' 



= sup ( Q 



2£ 



n 



for indices 0<l = n — k — 1 < n — 1. We recall that indeed Q is a symmetric matrix, so that 
More precisely, the error can be written 



EWu'' -V 



n\\2 



n-1 



k=o i,jes 

where for simplicity we use the notation := p(^^^) . . . P^^\ We compute for any i,j G 5, using 
the independence properties at different steps 

E((^fc)*(pW - g)Q2(n-l-fc)(p(fc) _ Q)A,),_^. 

= Yl lE[(A,)fc,,,(p('=) - Q)fc2,fei(Q'^""'"'^)fc2,fc3(^^'^ - Q)k„kAAk)k„j] 

ki,k2,k3,k4,&S 

= niAk)k,AAk)k„jmp^'^ - Q)k,Mip^'^ - Q)fc3,fc4](Q'^""'"'^)fc2,fc3- 

fci,A;2,fc3,fc4g5 

The observation is now that if k2 ^ k^, then the independence of the random variables for different 
nodes implies that 



(4.12) 



E[(P« - Q)fc,,fc,(p(*^) - Q),,3,fcJ = 0, 



since it is the covariance of two independent random variables - see (3.2). Moreover, when /c2 = ^3 
we see that {{Q^^~^~^^)*Q^"'~^~^^)k2,k3 only depends on n — A: — 1, due to invariance properties of 
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the equation. Therefore we rewrite the former expansion in the following way: 

n-l 



V 



ni\2 



k=0 ijes 



71-1 



(4.13) 



k=0 i,j&S 



k2&S 



n-l 

E 

k=0 



Q 



2{n~l-k) 



1,1 



E||(pW-Q)p('=-i)...p(0)u0||2,. 



k2,i 



We thus have to control ((5^^)i 



2l\ 



4.3 



gives a 



for any j ^ S. The following Lemma 
control of this expression. The first estimate means that the coefficients Q'j^ are approximations 
of the solution of the PDE at time 2i5t, at position j2, starting from the initial condition 0^^, with 
an error due to interpolation. The second estimate is fundamental in the proof of the Theorem, since 
it allows to introduce an additional factor 5x; however, we need to treat carefully the denominator. 



Lemma 4.3. There exists a constant C such that for any discretization parameters 5t 
5x = we have for any 1 < i < Mrp — 1 and for any < Ji, j2 ^ Ms — 1 



and 



(4.14) |Q,t,, - ^(t^ni^j, + B2m)\ < C—{1 + \ log{6t)\). 
Moreover, for any j £ S, we have for any 1 < i < Mt 

(4.15) ¥.<Pj{xj + B2m)<C- 



6t 
5x 



Remark 4.4. The singularities when — t- with a fixed 5x come from the use of regularization of 
the heat semi-group - when we consider the (f)j 's as initial conditions. 

For the second estimate (4.15), we make two important remarks. First, the constant C depends 

-\-oo: we have 

(j)j{x)dx = 5x^0. 



on the final time T, and we cannot directly let i tend to +oo; we have 

rXj+l/2 



lim ¥.(t)j{xj + B2m) 



Xj-l/2 



Second, from (4.15) we get for any ^ > and for any fixed 6t 

lim E(j)j{xj + B2est) = 0, 

Sx^O 

while we know that for a fixed 6x > and a fixed i, we have 

lim E(f)j{xj + B2m) = 4>j{xj) 
St^o 



1. 



These two behaviors are different and from (4.15) we see the kind of relations that the parameters 
5x and 5t must satisfy for obtaining one convergence or the other. □ 

Proof of Lemma 14.31 For any < i < 2Mt, we define 

Me = sup \iQ%j - E(j)j{xi + Bm)\, 
i,jes 

where {Bt)t>o is a standard Brownian Motion. 

We have Mq = 0, and by definition of Q we also have Mi = 0. 

We define some auxiliary functions Wj, for any index j £ S: for any a; S M and any t > 

Wj{t,x) =E(l)j{x + Bt). 
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Wj is solution of the heat equation, with periodic boundary conditions and initial condition For 
any t > 0, Wj{t, .) is therefore a smooth function - thanks to regularization properties of the heat 
semi-group - and since is bounded by 1 we easily see that we have the following estimates, for 
some constant C: 



(4.16) 



da^Wjit, .)||oo < and \\dl^W,{t, .)||oo < ^ 



We now prove the following estimate on the sequence (Af^): for any 1 < i < Mt — 1 

6x^ 



(4.17) 



M^+i < Mi + C 



i6t' 



The error comes from the interpolation procedure which is made at each time step. 
For any i,j G S, Markov property implies that 



{Q 



n,3 



kes 



= YQi,kiQ%d - EIoViWjii6t, .))(x, + Bst 

+ E[ior{Wj{e6t, .)) - Wj{m, .)]{xi + Bst). 

For the first term, we remark that it is bounded by M^; indeed we see that 

EI oV{Wj{£6t,.)){x^ + Bst) = YWjie5t,Xk)E(j)kixi + Bst) 



kes 

^<5i,fcE(/)j(xfe + B, 



est) 



kes 



To conclude, it remains to use the stochasticity of the matrix Q: entries are positive, and their sum 
over each line is equal to 1. 

The second term is bounded using the following argument: 



\\IoV(Wj{i5t,.))-Wjii6t,.)]\\oo < C5x'\\d'^^Wjii5t,.)\\oo < C 



.6x^ 

e5t' 



according to well-known interpolation estimates and to (4.16). 



From (4.17), using Mi = we obtain for any 1 < i < Mt 



Mt-1 2 
^^'^^^ E ^<t^^(|log(T)| + |log(<5t)|). 

k=l 



which gives the result, with a constant depending on T. 

Now we prove the second estimate of the Lemma. Thanks to the relation 0^,_|_i(a;) = (pki^ ~ '^2;) 
we see that the left-hand side does not depend on j E S; moreover we expand the calculation of the 
expectation using the periodicity of the function and relation definition of cj), the description of 
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its support as xj + IJfcezI^ ~ ^-^^ ^ + '^•^l- for 1 ^ ^ ^ Mt 



Ed 



1 



< 



1 



< 



< 



< C 



V2TTi6t 

1 

1 

5x 



k+Sx 

t 

— 5x 
k+Sx 



6x 



E 

kez 

E 

Y^C6xe 



—5x 
k+Sx 

I 

k—5x 



-k^/{2T) 



□ 



The estimate of Lemma 4.3 is now used in (4.13), and we obtain: 

n-l 



(4.18) 



E(<5 

k=0 



^2{ri-l-fe) 




|^(l + |.og(«)|)+^ 



(5x 



n 



(l + |log(50l) 



(5x 



+ 1 



1 - 



+ 1 



To conclude one more argument is necessary: we need to apply Proposition |4.1| in order to sum 
the variances. However this involves the quantity E|p(^~^) . . . P^^^u'^\'^i, which is badly controlled 
according to Proposition |3.4| for example, when 6t = 6x the accumulation only implies that 

6t 



E|p('=-i)...pW'u|2i <{l + C- 



CT 



u\ 



Ndx"^' 

for any k < j^. We recall that this bad behavior of the matrices p(") with respect to the h^-semi 
norm is a consequence of the independence of the random variables for different nodes, whereas this 
independence property is essential to get the improved estimate, since it allows to use the second 
estimate of Lemma 14.31 



Remark 4.5. Instead of considering Gaussian random variables which are independent with respect 
to the spatial index j, we could more generally introduce - like m j7] - a correlation matrix K, 
and try to minimize the variance with respect to the choice of K. Here we have chosen K as the 
identity matrix, so that the noise is white in space; the error hound (2.7) we obtain is a nontrivial 
consequence of an averaging effect due to this choice - see (4.12). A natural question - which is not 
answered here - would he to analyze the situation for general K: do we still improve the variance, 
and can we get more regular solutions? 

The solution we propose relies on the following idea: if above we could replace P^^-i) . . . p(0) 
with , we could easily conclude. Another error term appears, which is controlled by 1/A^ instead 
of l/y/~N. More precisely, independence properties yield for k>l 



(4.19) 



EIKP^'^) -Q)p(^-^)...p(0)n0||22 =E||(PW -Q)g'^n0|||2 



i/0||2 



+ E||(P('=) - Q) (^pC^-i) . . . p(o) - Q^'^ 
The roles of the differ ent terms are as follows. O n the one hand, the first term gives the part of size 



St 



thanks to Lemma 



4.3 



according to Corollary 4.2 and to Proposition 3.3 we have for any k>l 
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with kSt < T 
(4.20) EIKP^'^) 



On the other hand, the second term is now used to improve recursively the error estimate, since 
we have 



(4.21) E||(PW - Q) ( ^(^=-1) . . . - 



The independence of reahzations at step k gives the factor we remark that we cannot use the 
estimation of the one-step variance given by Corollary |4.2| otherwise we would need to control 
E|| (P('=-1) . . . - Q^) 



Using also (|4.18|) and (|4.19|) into (|4.13|), we see that 
(4.22) sup E\\u'' - v^l 



•""2 <C6tA^^ 



V 11^2. 



neN,n(5t<T " ^ 

The proof of the Theorem now reduces to the study of the following recursive inequalities, for p > 



with an initialization E^^^ 



C^, according to (4.11), with the notation B := {1 + 



St 



,0|2 



We 



remark that the control of the matrices P^^'^ and Q with respect to the /^-norm leads to another 
possibility for the initialization: E^^"^ = 2||n'^||^2; we observe that the recursion then yields the same 
kind of estimate. 

We finally easily prove that for any p > there exists a constant Cp > 1 such that 



(4.23) 



sup EWu^' 

nm,n5t<T 



V \\f2 2i 



(iVP 



AB 



and the proof of Theorem 2.2 is finished. 

o f:j2 

Remark 4.6. If we consider the equation ^ = with a viscosity parameter u > 0, the quantities 

A and B appearing in the proof are transformed into 



A], 



[1 + 



6x 



where the constant C does not depend on v. 



+ — (l + |log(5t)|)) and B, = {v^—)\uX^ 



5t 



The first change in the proof concerns the analysis of the one-step variance: in (|4.4|), the right- 
hand side is replaced by C{u5t + 5x^) 
same. 



We observe that the error due to interpolation remains the 



The second change concerns Lemma 4-3 where we use some regularization properties thanks to 
gaussian noise: when v goes to the estimates degenerates. 

As a consequence, we may observe that the estimate (2.7) gives a bound valid for a fixed value ofv, 
while (4.11) becomes more interesting when v is small compared with the discretization parameters. 



4.4. Accumulation of the interpolation error. To obtain Theorem |2.2[ it remains to control 
the deterministic part of the error of the scheme, without the discretization of the expectation with 
the Monte-Carlo method. We thus need to prove (2.6): 

for any n G N such that nbt < T, and for any j G N with < Xj = j6x < 1, we have 

(4.24) \u{n6t,Xj) 



< sup \uq{x) 



St 



xe[o,i] 



where u is the exact solution and where v"' is defined by (2.5). 
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Since \\u{n6t,x,) — f"||^2 < sup^ \u{n5t, xj) — v'j\, we easily obtain an estimate in the / -norm. 



Therefore, the conditions imposed on 6x and 6t by (2.7) are not restrictive, and can be seen as 
consequences of the semi-lagrangian framework. 



The proof of (4.24) in our context is as follows: using the exact representation formula and its 



discrete counterpart (2.5), we have 
u 



((n + l)6t, Xj) - = Eu{n6t, xj + Bst) - E vl(t)k{xj + Bs 

= ^{u{n5t,Xk) - v^)E(t>k{xj + Bst) 
kes 

+ E u{n5t, Xj + Bst) — u{n5t, Xk)4'k{'- 
\ kes 



where Bst is a Brownian Motion at time 6t. 
It is easy to see that 

I ^(n(n5t,Xfc) - v]^)Ecl)k{xj + Bst)\ < sup\u{n6t,Xk) - v]^\, 
kes 

and we see that the other term depends on the interpolation error: 

\E[u(ndt,Xj + Bst) — ''^^u{n6t,Xk)4>k{xj + Bst)]\ < sup \u{n6t,x) —I oVu{n5t, .)(x)\ 
kes ^e[o,i] 

<C5x^ sup \—-^{n5t,x)\<C5x^ sup |— — 2(0,a;)|. 
To conclude, we remark that for n = we have u(0, Xj) = v^. 

5. Numerical results and extensions 



5.1. Illustration of Theorem 2.2 The first numerical example we consider is a simulation of 
the solution of the heat equation in the spatial domain (0, 1) in periodic setting. We introduce the 
viscosity parameter so that the problem is 

(5.1) duit,x) ^ ^ 9M,x) ^ fo^. t>0,x£{0, 1), u{0,x) = uo(x)foT xe{0,l), 

ot ox^ 

with the boundary condition n(t, 1) = u{t, 0) for t > 0. For the numerical simulation of Figure [Tj we 
choose V = 0.01, and uq{x) = sin(27ra;). The exact solution satisfies x) = exp(— 47r'^i/t) sin(27rx). 
The discretization parameters are 6t = 5x = 0.01 and N = 100. 

The bound of Theorem 2.2 is illustrated with Figure [2| where we represent the error in logarithmic 



scales for different values of the parameters. 

We study the convergence of the scheme, with a numerical simulation which confirms the order 
of convergence with respect to the parameters 5t = 5x of the Monte-Carlo error. The final time 
is r = 0.1, the viscosity is = 0.1 and the initial condition is uq{x) = cos(27ra;). We compare 
the numerical solution vP with the exact solution; we only observe the Monte-Carlo error, which 



is dominant with respect to the deterministic part of the error according to Theorem 2.2 The 
mean-square error in the £^ norm is estimated with a sample of size 20. 

The error in Figure [2] is represented in logarithmic scales. The parameters 5t and 5x are equal 
and satisfy 5t = 5x = ^for the following values n = 50, 100, 200, 400, 800, 1600, 3200. Each line is 
obtained when we draw the logarithm of the Error as a function of \ogiQ[n)^ for a fixed value of 
N G {10,20,40,80}. The dot-line represents a straight-line with slope —1/2. 
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Figure 1. Solution at time T = 0.1 with 5t = 5x = l/N = 0.01 




1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 



Figure 2. Error for periodic boundary conditions when 5t = 5x = 1/n, in logarith- 
mic scales. 



This experiment confirms that the Monte-Carlo error is of order 1/2 with respect to the pa- 
rameters when 5t = 5x, as (2.7) claims. Indeed, the shift between the lines when N varies also 



corresponds to the size l/\/]Vof the Monte-Carlo error. 

5.2. The method for Dirichlet boundary conditions. We would like now show how it is 
possible to adapt our method in the case of Dirichlet boundary conditions. Let us consider the 

{0,1}. The 
= x + ^/vBi 

starting from the different points of the domain: If we define = inf {t > 0; G D"^}, then the 
solution satisfies 



equation (5.1), but with boundary conditions u{t,x) = for t > and x G dD 



representation formula then involves the family of the first-exit times of the process 



(5.2) u{t,x) =E[uo{X^)lt<r-]; 

the stochastic process is killed when it reaches the boundary. Note that this formula extends to 



more general PDF of the form (1.1) with the associated process (1.2). 



The numerical approximation becomes more complicated, since we also need an accurate approx- 
imation of the stopping times. This problem is well-known, and solutions have been proposed in |6j 
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and [9] for the computation of (5.2) at a given point x using time discretizations of the stochastic 
process Xf. 

In our case, we take advantage of the semi-lagrangian context to do a refinement near the bound- 
ary: for a discretization between the times t„ and tn+i, we introduce a decomposition of the domain 
into an "interior" zone and a "boundary" zone, with different treatments. In the boundary zone, 
we refine in time and use a subdivision of [n6t, (n + l)dt] of mesh size t < 5t and we use a possi- 
bly different value Nb for the number of Monte-Carlo realizations. Moreover, following [6j and [9], 
we introduce an exit test in the boundary zone, based on the knowledge of the law of exit of the 
diffusion process. 

In the interior part, less care is necessary and we can take t = 5t and Ni < for the size of the 
sample. 

We give in Figure [3] the result of investigations on the convergence of the method when Dirichlet 
boundary conditions are applied. We draw in logarithmic scales the error in terms oi n = l/6t = 
l/6x, with n = 50,100,200,400,800, with different values of the Monte-Carlo parameter Ni = 
10,20,40,80. We have chosen on the interval (—1,1) the initial function uq{x) = sin(7r^yi), with 
the viscosity ly = 0.1. The boundary zone is made of the intervals (—1, —0.9) and (0.9, 1), where we 
take r = 5t/10 and A^^ = lOiVj. The solutions are computed until time T = 0.1. Like in the case 
of periodic boundary conditions, the statistical error is dominant with respect to the other error 
terms; we compare with the exact solution, and to estimate the variance we use a sample of size 
100. 

The observation of Figure |3] shows that the Monte-Carlo error depends on the parameter 5t = 6x; 
the comparison with the "theoretical" line with slope —1/2 indicates a conjecture that the error is 
also of order 1/2, like for the periodic case. The shift between the curves for different values of 
corresponds in the error to a factor 1/^/N. 



-1.4r 




_2gl 1 1 1 1 1 1 1 

1.6 1.8 2 2.2 2.4 2.6 2.8 3 

log,„(n) 



Figure 3. Error for Dirichlet boundary conditions when 6t = 6x = 1/n, in loga- 
rithmic scales. 



5.3. The method for some non-Unear PDEs. We present a simple method to obtain approxi- 
mations of the solution of the viscous Burgers equation in dimension d = 2 

— + (n.V)u = uAu + f. 
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It is defined on the domain (—1, 1)^, with homogeneous Dirichlet boundary conditions - periodic 
ones would also have been possible. Compared with the situations described so far, we add a forcing 
term /, which may depend on time t, position x and the solution u. 

As explained in the Introduction, we construct approximations u" of the solution at discrete times 
n5t, introducing functions such that for any n > with the following semi-implicit scheme: 

(5.3) + (M'^.V)f"+i = i/At>"+i + r, 

for any time n6t < t < {n + l)6t and x ^ D. The initial condition is v^^^{n6t, .) = = v^{n6t, .). 
The discrete-time approximation then satisfies = uq and u"' = v'^{nSt, .). The forcing term here 
satisfies f^{t,x) = f{n6t,x,u'^{x)). 

On each subinterval [n6t, (n + l)St], we have 



\t,x) = E[v''+\n5t,Xf)lt<r- + / riX^)ds 

J nSt 

where the diffusion process X satisfies 

dXf = -u'^{Xf)dt + V2^dBt, X^st = ^■ 



The stopping times represent the first exit time of the process in the time interval [nSt, {n + l)5t]. 
Since v"'~^^{n6t, .) = u", the scheme only requires the knowledge of the approximations u". 

For the numerical simulations, we take the initial condition to be 0, and the forcing is f{t, x) = 
(— sin(7rf) sin(7ra:;) sin(7ry)^, — sin(7rt) sin(7rx)^ sin(7ry)). The viscosity parameter is = 0.001. The 
time step satisfies 5t = 0.02, and the spatial mesh size is 6x = 0.04. The "interior" zone is 
(-0.8, -F0.8)2, where Ni = 10; on the "boundary" zone, we have Ni, = 100 = lONi, and r = 0.002 = 
6t/lO. 

Both components of the velocity field u are represented in Figures |4] and [5] below at different 
times t = 0.5, 1, 1.5, 2. 
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Figure 4. Solution of the 2D Burgers equation at different times - first component 
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(a) U2,t = 0.5 



(b) U2,t = 1 





(c) U2,t = 1.5 (d) U2,t = 2 

Figure 5. Solution of the 2D Burgers equation at different times - second component 
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